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ABSTRACT 


This study examines mixing characteristics of double-diffusive convection for a wide 
range of fluids. Our approach involves Direct Numerical Simulation (DNS) utilizing de- 
aliased pseudo-spectral method. To expedite these simulations the numerical algorithm 
was parallelized using Message Passing Interface (MPI) calculations in both two and 
three dimensions. A theoretical model of equilibrium double-diffusive transport is 
presented, which emphasizes the role of secondary instabilities of salt fingers in 
saturation of their linear growth. Theory assumes that the fully developed equilibrium 
state is characterized by the comparable growth rates of primary and secondary 
instabilities. This assumption makes it possible to fonnulate a simple and efficient 
algorithm for computing diffusivities of heat and salt as a function of the background 
property gradients and molecular parameters. The model predicts that the double- 
diffusive transport of heat and salt rapidly intensifies with decreasing density ratio. 
Fluxes are less sensitive to molecular characteristics, mildly increasing with Prandtl 
number (Pr) and decreasing with diffusivity ratio (x). Theory is successfully tested by a 
series of direct numerical simulations which span a wide range of Pr and x. 

Double diffusion occurs on the micro-scale and computer technology is just now 
reaching the processing speeds needed to fully resolve this complex phenomenon in three 
dimensions. In addition to the well-known “salt fingering” within the oceans, double 
diffusion occurs in many other statically stable regions, both terrestrially and beyond our 
atmosphere. Understanding this phenomenon could prove essential as we continue to 
discover new worlds and new areas within our galaxy, including here on terra finna. 
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I. INTRODUCTION 


Double diffusion is the instability of a stratified fluid at rest whose density is 
determined by two components diffusing at different rates. Such a configuration can be 
unstable even if the density of the fluid is increasing downwards. The resulting double- 
diffusive convection has long been recognized as a significant, and in many cases 
dominant, mixing process in the ocean, where wann and salty water is often located 
above cold and fresh. In this case, the faster diffuser (temperature T) is stabilizing and the 
slower diffuser (salinity S) is destabilizing, resulting in the salt finger form of double- 
diffusive convection (Figure 1), which is the main focus of this discussion. 
Equation Section (Next) 

Salt fingers are mere millimeters in diameter and can grow to a length of over 20 
centimeters (Pickard 1990). Double-diffusive convection interactions were originally 
noted by W. S. Jevons (Jevons 1857) and were the first recorded occurrences of this 
phenomenon. At the time, Jevons was trying to model cirrus clouds when he 
inadvertently discovered double-diffusive convection within the laboratory environment. 
In his entry to Edinburgh and Dublin Philosophical Magazine and Journal of Science in 
1857, Jevons used a warm sugar solution over a fresh water solution in a lab experiment 
and reported below: 

The parts of these strata, however, which are immediately in contact, soon 
communicate their heat and tend to assume a mean temperature; and it is 
evident that whenever this is the case, the portions of liquid containing 
sugar must always be slightly denser than those that are pure, and must 
consequently sink below and displace the latter. 

We shall thus have portions of the upper stratum continually sinking into 
the lower, and corresponding portions of the lower rising through the 
upper; and this movement, as the experiment demonstrates, takes place by 
an interfiltration of minute, thread-like streams. 

The first physical model of salt fingers was developed by Dr. Melvin E. Stern 
(Stern 1960) more than a century later, which marks the beginning of the modern theory 
of double-diffusive convection. This phenomenon is better understood today due to the 

use of high performance computers (Radko and Stern 1999; Kimura and Smyth 2007; 
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Prikasky 2007; Caro 2009) important insights were brought by laboratory studies 
(Krishnamurti 2003) and field measurements with in the ocean (Gargett 1982). While it is 
still not fully understood how these minute phenomena impact larger scale circulations, it 
is widely believed that they may hold the key to understanding the driving forces behind 
the general circulation of the ocean and therefore, the Earth’s climate. 

In addition to oceanographic applications, compositionally driven double- 
diffusive convection affects heat and material transport in a variety of other geo- and 
astrophysical fluid systems, from magmatic melts (Tait and Jaupart 1989) to the interiors 
of giant planets and stars (Guillot 1999; Vauclair 2004; Charbonnel and Zahn 2007; 
Stancliffe et al.„ 2007). 

A. SALT-FINGERS IN THE OCEAN 

Salt fingers occur in nature in statically stable environments (Figures 2 and 3), 
where temperature and salinity perturbations diffuse at different rates causing fingers to 
form within a body of water. Salt fingers are common in the main thermocline, where 
water tends to be wanner and saltier near the sea-surface due to a combination of surface 
evaporation and heating. Individual parcels displaced downward lose their thermal 
properties faster than salinity, causing them to be denser than the surrounding fluid 
(Figure 1) and continue to sink. This causes the parcel to continue descending (Stern 
1960; Schmitt et al.„ 2005), which results in vigorous double-diffusive convection 
discussed in this study. 
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Salt Fingers 



Figure 1. Schematic diagram illustrating the physical mechanisms of double-diffusive 

convection. 

An example of this process can be found at the output of the Mediterranean 
through the Straits of Gibraltar. The outflow from the Mediterranean Sea is much saltier 
and wanner than the Atlantic Ocean water at depth (approximately 2400 meters) causing 
double-diffusive interactions deep underneath the surface of the ocean. 

It is also believed that salt fingers contribute to the ocean’s large scale mixing and 
could be important to climate models (Gargett and Schmitt 1982). Today’s climate 
models cannot resolve such minute scales within a reasonable forecast time, so finding 
physical parameterizations that can be fed into these models can greatly improve their 
accuracy. Understanding the impact of these processes on large scale circulations can 
also help to provide better ice melt forecasts, as well as higher resolution global and 
regional oceanographic forecast models across all spectrums within oceanography. 
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Figure 2. Diagram illustrating the geographic distribution of diffusive-convection, 
which occurs in approximately 15% of the world’s oceans. Unlike salt fingering, 
this phenomenon occurs when cold fresh water is located above warm and salty. 



Figure 3. Diagram illustrating geographic locations of salt fingers, which occurs in 

approximately 30% of the world’s oceans. 

B. BEYOND OCEANOGRAPHY 

By extending our knowledge of oceanographic double-diffusion and changing 
background parameters we can get a “feel” of how double diffusive convection changes 
in other environments. 

External planetary systems, such as the atmospheric makeup of planets within our 
solar system, are represented by a much lower Prandtl number (Pr) than in our oceans, 
(Pr = 10° verses Pr = 7). Estimates show that for Jupiter’s atmosphere (radiation driven), 
Pr « 10' 4 and for the metallic interior (conduction driven), Pr « 10' 2 (Bagenal et al.„ 
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2004). Decreasing Pr from the oceanographic value of 7, allows us to build a better 
understanding of how double diffusive convection behaves within these foreign bodies. 

The interior temperature and Helium within our own sun and other stars is also a 
natural area where double diffusive convection occurs as interior temperature and Helium 
diffuse at separate rates. Metallic lingers are theorized in these systems to consist of 
metallic “blobs” that descend similar to the oceanographic warmer denser fluid, with heat 
diffusing much quicker than the heavy metal composition within the fluid (Vauclair 
2003). Another manifestation of active salt lingering in nature, the high Prandtl number 
within magma flows, has close similarities to oceanographic salt lingers (Tait 1989). 

Sugar-salt interactions have been modeled in the laboratory environment starting 
with Jevons in 1857 with pure water and, hydrochloric acid and white sugar (Jevons 
1857). The sugar-salt solution is different to prior cases, in both the Prandtl number and 
the diffusivity ratio. However, we include this case so as to achieve a more complete 
understanding. 

While many regimes are currently beyond the modeling capabilities of even for 
the most powerful super computers, in this study we strive to gain a quantitative 
understanding within the accessible limits of Pr, x, and R p values. Many regimes we wish 
to explore are outside this reach and we are pushing the limits of what is possible with 
today’s computer processing speeds. 

C. EXTANT MODELS OF EQUILIBRIUM TRANSPORT 

A fundamental question in double-diffusive convection theory concerns the 
equilibrium transport of temperature, salinity, chemical tracers, and momentum. The 
quantification of double-diffusive fluxes and their dependencies on the background 
temperature and salinity gradients—so-called flux-gradient laws—is an essential step in 
linking the microstructure dynamics with larger scales of motion. This problem has 
motivated numerous laboratory (Lambert and Demenkow 1972; Griffiths and Ruddick 
1980) and field (Schmitt et al.„ 1987; Schmitt et al.„ 2005) experiments, as well as 
theoretical (Stern 1969; Kunze 1987) and numerical (Shen 1995; Stern et al.„ 2001; 
Kimura and Smyth 2007, 2011; Traxler et al.„ 2011; Stellmach et al.„ 2011) models. 
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While the linear stability theory of double diffusion (Stem 1960; Baines and Gill 1969) is 
well developed and fully understood, it does not explain the saturation of primary double- 
diffusive instabilities—the problem of equilibrium double-diffusive transport is 
complicated and principally nonlinear. 

Numerous attempts have been made to deduce flux-gradient laws from first 
principles. The first and perhaps the most influential idea was introduced by Stern 
(1969), who suggested that the linear growth of salt fingers is arrested when the Stern 
number 

£ — ~ PFsi\m 

~ v(aT - [3S z ) 

reaches 0(1). F T dim and Fsj, m here are the dimensional temperature and salinity fluxes; T, 
and S_ are the vertical gradients; a, fi are the thennal expansion and haline contraction 
coefficients respectively; v is the molecular viscosity. Stern’s physical argument is 
compelling. When A exceeds unity, the unbounded salt finger system becomes unstable 
with respect to collective instability, a term used to describe the spontaneous excitation of 
gravity waves by salt fingers. Stem suggested that collective instability could dismpt salt 
fingers, thereby arresting the development of primary instabilities. 

At first glance, various pieces of indirect evidence seemed to validate Stem’s 
hypothesis. Oceanographic measurements (Hebert 1988; Innoue et al.„ 2008) tend to 
produce 0(1) values of A. Kunze (1987) gave an alternative argument in support of the 
Stern number constraint. He showed that A = 1 is equivalent to the Richardson number 
(Ri) of % based on scales of individual salt fingers. Kunze suggested that the well-known 
criterion for the instability of horizontal, inviscid and nondiffusive parallel flows 
(Richardson 1920; Howard 1961; Miles 1961), i.e. Ri < %, can be extended to viscous, 
diffusive and vertically oriented fingers. 

However, closer inspection has revealed certain inconsistencies in the Stem- 
Kunze hypothesis. For instance, in many laboratory experiments, salt and sugar replace 
heat and salt as buoyancy components (Stern and Turner 1969; Lambert and Demenkow 
1972; Griffiths and Ruddick 1980; Taylor and Veronis 1996; Krishnamurti 2003). 
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Typical values of Stern number in this case are extremely low. Lambert and Demenkow 
(1972) report Stem numbers as low as A = 2-10' , casting some doubt on the generality of 
the Stern number constraint. On the theoretical side, a critical advance was made by 
Holyer (1984), who performed a formal linear stability analysis for vertical steady salt 
fingers. She discovered direct, relatively small-scale — comparable to the salt linger 
width — secondary instabilities, which typically grow much faster than the collective 
instability modes. Unlike collective instabilities, Holyer modes grow regardless of the 
(finite) amplitude of salt fingers and thus regardless of the specific values of A. 
Numerical simulations (Shen 1995; Radko and Stern 1999) confirmed Holyer’s ideas by 
demonstrating that the equilibration of salt fingers occurs by means of the nonlinear 
interaction of primary instabilities with Holyer modes. Furthermore, the eddy 
diffusivities of heat and salt from simulations (Stern et al.„ 2001) and available 
observations (St. Laurent and Schmitt 1999) monotonically decrease with the density 
ratio, whereas the opposite trend was expected based on the Stern-Kunze constraint. The 
most recent and well-resolved simulations of the heat-salt system (Traxler et al.„ 2011) 
indicate that as the density ratio increases from 1.2 to 10, the Stern number in the 
equilibrium flow reduces by more than two orders of magnitude, from A = 76 to ^4 = 0.6. 

Concerns with regard to the Stern number constraint have motivated the 
development of alternative models for the equilibrium transport. While the general 
analytical description of double-diffusive transport is still lacking, several cases have 
already been successfully treated. Weakly nonlinear instability theory describes salt 
finger patterns, dynamics, and transport characteristics for the parameters near the point 
of marginal instability (Proctor and Holyer 1986; Radko and Stern 1999, 2000; Stern and 
Simeonov 2004; Radko 2010). Promising attempts to formulate physically based 
parameterizations of salt finger transport include the mixing length model (Stem and 
Simeonov 2005), a double-diffusive version of the upper bound theory (Balmforth et al., 
2006), a second-order closure approach (Canuto et al., 2008), and various 
phenomenological models (Shen 1995; Merryfield and Grinder 1999; Radko 2008). We 
propose what appears to be, to date, the most general algorithm for estimating the 
equilibrium salt finger transport. The model is based on the properties of secondary salt 
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finger instabilities and can be applied to a variety of parameter regimes, including small, 
moderate, and large Prandtl numbers, and a wide range of diffusivity ratios. 


D. OUR APPROACH 

Our story is a new variation on an old theme—a theme that appears, explicitly or 
implicitly, in almost all theories of equilibrium double-diffusive convection (Stern 1969; 
Holyer 1984; Kunze 1987; Stern and Simeonov 2005; Radko 2010). We revisit the idea 
of a competition between linear mechanisms involved in the growth of salt fingers and 
the disruptive action of their secondary instabilities. At the same time, we strive to avoid 
some debatable assumptions of earlier models, such as the link between secondary 
instabilities and Stern number. The competition between primary and secondary 
instabilities is concisely phrased in the growth rate balance assumption: 

( 1 ) 

where X\ is the typical growth rate of linear salt fingers, X 2 is the growth rate of their 
secondary instabilities, and C is a dimensionless order one quantity that can be calibrated 
on the basis of simulations or experiments. The primary growth rate k| is detennined by 
the background gradients of temperature and salinity and by molecular characteristics 
(diffusivities and viscosity). The secondary instability X 2 is also affected by these 
quantities, but, additionally, it depends very strongly on the amplitude of primary salt 
fingers. Thus, for any given background parameters and a value of C, the growth rate 
balance Equation (1) implicitly determines the equilibrium amplitude of salt fingers. 

The physical reasoning behind Equation (1) is straightforward. As initially weak 
salt fingers, emerging from random small scale perturbations, grow in time, they start to 
develop secondary (Holyer 1984) instabilities. Unlike the primary ones, the growth rate 
of secondary instabilities monotonically increases with the amplitude of fingers. At first, 
the growth of secondary instabilities is too slow to inflict any significant damage on 
growing primary modes—the evolution of small amplitude fingers is adequately captured 
by the linear theory. However, at some point, the growth rate of secondary instabilities 
exceeds the primary growth rates. As a result, the secondary instabilities start to gain in 
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magnitude, rapidly reaching the level of primary modes. Their interaction nonlinearly 
suppresses the growth of salt lingers. At this stage, the system reaches statistical 
equilibrium. 

In this study, condition Equation (1) is used as the basis for a simple algorithm to 
determine the equilibrium amplitude of salt fingers. The growth of primary modes (ki) is 
evaluated from the linear instability theory (Stern 1960; Baines and Gill 1969). For any 
given amplitude of salt fingers, the growth rate of secondary instabilities (X 2 ) can be 
computed using elements of the Floquet theory as in Holyer (1984). The amplitude of 
fingers is iteratively adjusted until X\ and Xj satisfies the growth rate balance Equation (1) 

. The model predictions are then conveniently expressed in terms of the equilibrium eddy 
fluxes of heat and salt. 

The use of DNS has been essential to the success of our research. DNS allows us 
to (i) calibrate our theoretical model, (ii) validate our reference base on the growth rate 
balance Equation (1), and (iii) expand the sensitivity of our solutions to the background 
parameters. From crude calculations of man-hours spent on this project it would take 
91+ years (over 800,000 hrs) to run these calculations on a single computer processor 
continuously. Many of these experiments are run on multiple parallel processors, 
sometimes 128 working simultaneously and still taking longer than a week to complete. 
As computing power increases, this processing time will decrease and more robust 
calculations and data can be collected. 

This thesis is organized as follows. In Chapter II, we present DNS, focusing our 
inquiry on the equilibration of primary instability. A theoretical model of equilibration is 
presented in Chapter III. We compute the temperature and salinity fluxes as a function of 
density ratio for the oceanographic case (Pr=7, r=0.01) and successfully test (Chapter IV) 
the theoretical prediction by DNS. In Chapter V, we also explore a broader parameter 
range, including the low Prandtl number regime, relevant for astrophysical applications, 
and high Prandtl numbers, relevant for magmatic melts. Each case is compared with the 
corresponding DNS. We summarize and draw conclusions in Chapters VI while 
discussing future possibilities in Chapter VII. 
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II. PRELIMINARY CALCULATIONS 


The simplest and most direct approach to the analysis of equilibrium of salt 
fingering involves Direct Numerical Simulation (DNS). We used a de-aliased pseudo- 
spectral method described and first used for two-dimensional (Stern and Radko 1998) and 
three dimensional salt-fingering simulations (Radko and Stern 1999). In order to 
expedite the simulations the numerical algorithm was parallelized using a Message 
Passing Interface (MPI) algorithm in both two and three dimensions. Two dimensional 
model runs were performed on Anastasia, a high-performance computer cluster at Naval 
Postgraduate School (NPS). Three dimensional model runs were perfonned on the 
shared high-performance computer cluster Hamming, also at NPS, and Einstein, a U.S. 
Navy Cray XT5 high-perfonnance computer cluster. The resources utilized, varied from 
4 processors on Anastasia up to 128 simultaneous processors on Einstein for the most 
highly resolved cases in heat salt throughout R p values of 1.1, 1.3, ..., 2.9 (Tables 2 
through 6). We resolve a total of ten fingers; 5 warmer, more dense and 5 cooler, less 
dense. However, in order to resolve the salt finger scale, these runs required anywhere 
from 3 days to 4 weeks to complete. With increased computing power, this time 
limitation to process data will decrease in future years. In turn, this will allow future 
graduate students the flexibility to run more data sets over longer periods and thereby 
address even more challenging problems in double-diffusive convection. 


Following Radko and Stern (1999), we separate the temperature and salinity fields 
into the basic state , representing a uniform vertical gradient, and a departure 

(T, S) from it. The two-dimensional Boussinesq equations of motion are expressed in 


terms of T,S and nondimensionalized using / = 


f k T v V k 


—, —, and ^° V ^ T as the 


/ k T 


l 1 


\S aT :J 

scales of length, velocity, time and pressure respectively. Here, ( kr,ks ) denote the 
molecular diffusivities of heat and salt and po is the reference density used in the 
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Boussinesq approximation. The expansion/contraction coefficients are incorporated in 
(T,S) and aTzl is used as the scale for both temperature and salinity perturbations, 
resulting in 


dT 



+ V- 

dt 


8S 



+ V- 

dt 


1 | 

f— 


w 


R 


Pr 


A 


— v + v • Vv 
dt 


= —Vp + (T- S)k + Vv 


J 


V-v = 0, 


( 2 ) 


where v = (w,v,w) is the velocity vector. This system is unstable with respect to salt 
lingering (Stem 1960) for 

\<r p <\/t. (3) 

The key nondimensional numbers governing the evolution of system Equation (2) 
are the Prandtl number p r = — , the diffusivity ratio r = —, and the background density 

kj k T 

aT 

ratio R = — ==-. We also assume that, in the absence of large-scale structures, fluxes are 

p ps z 

independent of the nondimensional parameters related to the domain size (e.g., the 
Rayleigh number). The local flux-gradient laws are commonly used to parameterize the 
effects of salt fingering on the oceanic circulation and successful attempts have been 
made to validate the concept of an unbounded T-S gradient in numerical simulations 
(Radko and Stem 1999). 

To gain a preliminary understanding of the mechanics of equilibration (R p ,Pr,t), 

Equation (2) was solved numerically. We assume triply periodic boundary conditions for 

T, S, p, and v and integrate the governing equations using a de-aliased pseudo-spectral 

method described in Radko and Stern (1999). In the following calculation, we use 

parameters relevant in the oceanographic (heat/salt) context: z=0.01 and Pr=7. The 

overall density ratio is R p =l .9, which is representative of the main Atlantic thermocline. 

The size of the computational domain (L x = 38.2, L y = 38.2, L z = 76.3) is equivalent to 
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(5x5x10) linearly fastest growing finger wavelengths ( d ). The flow was resolved by a 
uniform mesh with (N x x N y x N z ) = (256 x 256 x 512) elements, and the model was 
initialized from rest by a small-amplitude random computer-generated initial ( T,S ) 
distribution. After a few characteristic growth periods, active statistically steady double- 
diffusive convection was established. Figure 4a shows a typical instantaneous (t= 20) 
temperature field in the initial stage of linear growth. 



Figure 4. Equilibration of salt fingers in the numerical experiment with Pr=7, x = 0.01, 
R p = 1.9. Three-dimensional instantaneous temperature fields are shown for a) 
the early stage of linear growth at t =20 and b) the fully equilibrated state at t= 50. 
Red/green color corresponds to high values of T and low values are shown in blue 
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As expected from the linear stability theory, the most rapidly growing 
perturbations take the form of the vertically uniform salt fingers known as elevator 
modes. Figure 4b presents the fully developed equilibrium regime (t=50). While the 
previously dominant elevator mode is still visible, it is now comparable in magnitude to 
the secondary instabilities, which act to twist the vertical fingers, adversely affecting their 
growth. This stage is characterized by the irregular transient patterns in the temperature 
field—a dramatic consequence of the nonlinear interaction between the elevator modes 
and their secondary instabilities. 



Figure 5. Time record of the temperature (solid curve) and salinity (dashed curve) 
fluxes for the calculation in Fig 4. Equilibrium point is marked by red line. 


The transition to the statistically steady regime is illustrated in Figure 5, which 
shows the evolution of the spatially averaged downward temperature and salinity fluxes 
(F t ,F s ) = (-wT,-wS) - the sign convention is dictated by considerations of convenience. 
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Note that the heat flux in our system of nondimensionalization is F T = Nu -1, where Nu 
is the Nusselt number—the ratio of the total heat flux to the (much weaker) molecular 
diffusive flux. The initial stage of the simulation (t<30) is characterized by the 
exponential growth of fluxes, which is followed by their equilibration. After 
equilibration, the intensity of salt fingering remains statistically steady. 

The nondimensional fluxes ( F T , F s ) can be converted to the dimensional eddy 
diffusivities of heat and salt (K Tdim ,K Sdim ) using 

^TUm = F T k T , K Siim = F s k T R p . (4) 

The time mean T-S fluxes, averaged over the second half of the experiment in 
Figure 5, are F T =40.88 , F s =75.88, which corresponds to the dimensional diffusivities 
of heat and salt of 

K T6im =5.6-10- 6 mV 1 , K Sdim = 2.0 • 10- 5 m V 1 (5) 

These values are broadly consistent with the observational estimates (St. Laurent 
and Schmitt, 1999) for the density ratio of R p = 1.9. 

In Figure 6, we plot typical horizontal (top panel) and vertical (lower panels) 
sections of temperature (left) and salinity (right) for the fully developed equilibrium 
regime. Temperature and salinity patterns are correlated and both reveal vertically 
elongated—but far from uniform—horizontally isotropic salt fingers. Salinity sections 
are characterized by very fine, relative to temperature, spatial scales, which is a direct 
consequence of the low diffusivity ratio. 
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Figure 6. Horizontal (upper panels) and vertical (lower panels) sections of temperature 
(left) and salinity (right) for a typical fully equilibrated state of the calculation in 

Fig 5. 
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III. THEORETICAL MODEL 


Guided by physical reasoning (Section D of Chapter I) and by direct numerical 
simulations (e.g., Chapter II), we now attempt to predict the equilibrium level of salt 
fingers from the growth rate balance Equation (1). The growth rate equation for the 
primary salt finger instabilities (Stem 1960) in our nondimensional units reduces to 

A 3 + [l + r + Pr]A : 2 A 2 + [(r + Pr+ rPr)/t 4 + Pr (l - R p ' )] A + r Pr k 6 -Pr k 2 (i?; 1 - r) = 0, ^ 

where k is the horizontal wavenumber of the vertically uniform elevator modes. For our 
purpose it is not necessary to consider inclined salt fingers since the largest growth rates 
are attained by the vertically oriented elevator modes (e.g., Baines and Gill 1969). For 
each (R p ,Pr,r), we compute the largest growth rate (A|) by maximizing the solutions of 
Equation (6) with respect to k as in Schmitt (1979, 1983). Thus, Ai is uniquely 
determined by the governing parameters: 

4 =/((R ,,Pr,r). (7) 

While the maximal growth rates of primary instabilities are identical in two (x,z) 
and three (x,y,z) dimensions, the secondary instabilities differ substantially. Therefore, 
the algorithms for computing An in two and three dimensions will be discussed separately. 


A. TWO-DIMENSIONAL FORMULATION 

In two dimensions, the governing equations are conveniently written in terms of 

the streamfunction 1 // , such that (u,w) = ( --—j, resulting in 

V dz dx J 


dT dw , 

— + J(y/,T)-\—— = V 2 T 

dt dx 

dS 1 dw 

— + J{y/,S) +-— = zWS 

dt R dx 


( 8 ) 


d ^ 0 

— V-i// + J(i//,'V~i//) = Pr 
dt 


A (r _ 5 ) + vy 

dx 


Our next step is to use Equation (8) to examine the stability of the basic state 
representing the elevator modes: 
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T = A t sin (kx) 

< S - A s sin(foc) (9) 

yJ - A w cos(Gc) 

It should be mentioned that the secondary stability problem can be formulated in 
different ways. Holyer (1984) has chosen the wavenumber k in Equation (9) to 
correspond to the solution with zero primary growth rate (k 0 ). This restriction defines a 
well-posed stability problem with the regular basic steady state. There is only one reason 
to question this formulation. It is apparent from DNS (Stern et al.„ 2001; Traxler et al.„ 
2011) and even indicated by some oceanographic field measurements (Gargett and 
Schmitt, 1982) that the horizontal wavenumber of dominant fully developed salt lingers 
is more accurately approximated by k max , the wavenumber corresponding to the 
maximum growth rate, than by k 0 . Use of k=k max in Equation (9); however, results in the 
time-dependent basic state, and therefore the conventional methods of stability analysis 
cannot be applied directly. This problem is not uncommon and not insurmountable. 
Notable examples of analyzing such flows come from studies based on the Kolmogorov 
model—the sinusoidal parallel flow in viscous fluid (Sivashinsky 1985; Manfroy and 
Young, 1999 2002; Balmforth and Young 2002, 2005). Since the unforced viscous 
parallel flow would inevitably decay, the Kolmogorov model circumvents this problem 
by introducing artificial forcing in the momentum equation that maintains the steady 
state. Variations on the same principle, often referred to as the quasi steady state 
approximation or a “frozen flow” method, have been successfully applied to numerous 
stability problems (e.g., Lick 1964; Robinson 1976) including salt fingers (Kimura and 
Smyth 2011). 

A choice has to be made at this point between the approach of a purist, insisting 
on formal mathematical consistency, and that of a practically-oriented researcher, more 
concerned by the consistency with patterns realized in nature. We have examined both 
models, k=ko and k=k max , and the results are qualitatively similar. The prediction of the 
k=k max model is closer to the corresponding DNS, and therefore we limit the following 
discussion to the configuration based on the fastest growing solution. We assume that in 
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the regime close to equilibrium the amplitude Equation (9) is maintained at a quasi¬ 
steady level for k=k max and separate the solution into the dominant basic state and small 
perturbation: 

{T,S,y/) = {T,S,w)+(T',S',y/'). (10) 


Equation (8) is linearized as follows: 


dT' , , . .ar , , /( .dur' dy/' - , 

- Ak sin(fcc)- A r k cos (he) +-= V~T 

dt V dz T dz dx 

dS' dS' . , ,, .dy/' 1 dy/' _ 2o , 

- A k sin (kx) - A^kcos(kx) -h-= iV~S 

dt v dz s dz R dx 


(ID 


— Vi// - A ¥ ksm(kx) — - A y/ k sin(kr)^ = Pr 


—(T'-s')+vy 

dx 


To examine the stability properties of the linear system Equation (11), we use the 
Floquet technique in which the perturbation is sought in the following form: 



N 

f T \ 

1 n 

S' 

= exp {ifkx + imz + At) ^ 

s. 

y j 

n=-N 



cxp(inkx ), 


( 12 ) 


where k is the growth rate, m is the vertical wavenumber, and/is the Floquet coefficient, 
which controls the fundamental wavelength in x. Substituting Equation (12) into the 
linear system Equation (11) and collecting the individual Fourier components allows us 
to express the governing equations in the matrix form: 




'Av| 

II 

(13) 

where 


d, — (T_ N ,S_ N ,i//_ N ,T_ N+i ,S_ N+ i,i//_ N+l , ... ,T n ,S n ,i// n ) , 

(14) 

and A 

is 

the square matrix whose elements are 

functions 


of(k,m,f,R p ,Vx,T,N,A T ,A s ,A v ). 

The growth rates of the normal modes of the linear system Equation (11) 
correspond to the eigenvalues of A . For each set of governing parameters, we determine 
the eigenvalue with the maximum real part, which represents the fastest growing mode. 
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The value of this growth rate is then maximized with respect to m and f Note that the 
symmetry and periodicity properties of our system are such that the Floquet coefficient/ 
needs to be varied only within the range 0 < f < 0.5. 


The coefficients (A t As,A v ) of the normal mode Equation (9) are connected, for 
k=k max , by the following relations: 


A T (A + k 2 ) = k A 

T \ 1 max J max y/ 


(15) 


and therefore the amplitude vector (A T ,As^4 v ) can be immediately determined from any 
component—for instance from the temperature amplitude (At). As described below, the 
solution rapidly converges with increasing resolution (N—> oo). Thus, the fastest growing 
secondary instability is determined by four parameters: 

A 7 =A 7 (R p ,Pr,r,A T ). (16) 


Our next step is to solve the growth rate balance Equation (1). For each set of 
governing parameters (R p ,Pr,i:), we determine the primary growth rate / using Equation 
(6). Suppose now that the value of C in Equation (1) is known. Then, by varying the 
amplitude of the normal mode Aj we can readjust Equation (16) until the growth rate 
balance is satisfied. The dependence A 2 (d T ) takes the form of a monotonically increasing, 
nearly linear relation. For Aj = 0, the secondary instability problem becomes identical to 
that of primary instability and therefore X 2 = h ■ A corollary of this observation is that 
solutions of Equation (1) exist only for Of as indicated by the schematic in Figure 7. 
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Figure 7. Illustration of the algorithm for predicting the equilibrium amplitude of salt 
lingers. The growth rate of the secondary instability (X 2 ) monotonically increases 
with the amplitude of salt fingers (A r ). We hypothesize that the equilibrium (A*) 

is reached when X 2 =C Xi. 


An iterative procedure for solving the growth rate balance Equation (1) was coded 
in MAPLE and it typically required 7-8 iterations for the model to converge, within a 
negligible relative error of 10' 5 , to the sought-after solution^ = T* . The vertical fluxes 
are then reconstructed using Equation (9) and Equation (15) as follows: 


f t = wt = '-4 (A+C.) 


2 f? 


(4 + r ^max ) 


(17) 
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Figure 8. Heat flux in the theoretical two-dimensional model for various values of C. 

In Figure 8, the predicted heat and salt (Pr=7, x=0.01) fluxes are plotted as a 
function of density ratio (R p ) for various values of C. Fluxes increase with C and rapidly 
reduce with R p . The calculations in Figure 8 were perfonned for N= 32. However, the 
specific number of Fourier harmonics involved in the Floquet calculation is not 
particularly important. The rapid convergence of the solution with increasing N is 
indicated in Table 1. The very crude truncation with N= 2 results in 19% error, whereas 
the N— 4 error is less than 2%. 
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N 

Ft 

F s 

F t —F t 

rrr - T T V=32 

UA A — . 

F t \ 

t \N= 32 

2 

29.55387 

50.23965 

0.1899 

4 

25.28658 

42.98554 

0.0181 

8 

24.83962 

42.22573 

0.0001 

16 

24.83668 

42.22074 

err < 1CT 5 

32 

24.83668 

42.22074 

NA 


Table 1. Convergence of the Floquet based algorithm in a two-dimensional 
calculation for (R p ,Pr,x) = (1.9,7,0.01). As the number of spectral harmonics ( N ) 
increases, the heat flux (Ft) rapidly approaches its limiting (N—»oo) value. 


B. THREE-DIMENSIONAL FORMULATION 

Numerical simulations (e.g., Stern et al.„ 2001) reveal significant differences, in 
magnitude and dynamics, between two- and three-dimensional salt fingers. Three- 
dimensional fluxes are systematically higher, by a factor of 2-3, than their 2D 
counterparts. In this section, we attempt to extend the growth rate balance theory to three 
dimensions. In addition to the natural interest in three-dimensional dynamics—after all, 
we live in a three-dimensional world—we are also motivated by the expectation that the 
results will help to rationalize the difference in the intensity of salt lingers in 2D and 3D. 
Since the following theory is analogous to the two-dimensional case (Section A of 
Chapter III), it is presented in abbreviated form. 

The essential difference between two- and three-dimensional dynamics stems 
from the more complicated structure of the salt finger elevator mode in 3D: 
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T = A t exp(/l/)^(jc, y) 

< S = A s exp(A 1 t)^(x, y) (18) 

w = A w Qy^{\t)(j){x,y), 


where the horizontal planform function (j> satisfies the Helmholtz equation V 2 H t/> + k 2 cj) - 0 . 
In 3D, identical growth rates can be attained by various horizontal planforms. For 
instance, even if we limit our analysis to rectangular planforms cf> = cos(k x x)cos(k v y ), all 


possible configurations with 


+ K = k, 


2 

max 


(19) 


yield the maximum growth rates (Ai) and therefore are a priori equivalent. The 
degeneracy of linear theory with respect to the planform selection is well known and can 
be resolved by invoking fundamentally nonlinear considerations. Proctor and Holyer 
(1986) and Radko and Stern (2000) examined the problem for a “bounded” configuration, 
in which finger layer was limited by rigid horizontal surfaces—the setup analogous to the 
classical Rayleigh convection problem. The former (latter) study suggested the tendency 
of the planform to evolve spontaneously to two-dimensional rolls (square cells). For the 
unbounded regime, some guidance regarding the planform selection was provided by 
diagnostics of numerical simulations in Radko and Stern (1999), which indicate the 
preferred planform is neither square nor a roll, but, rather, some combination thereof, 
perhaps best described by the planfonn function of the following type: 

^ = cos(£ max x) + -^-cos(& max y), (20) 

which was used in all calculations described below. It should be mentioned, however, 
that for our purpose—determination of the equilibrium transport characteristics—the 
planform choice turns out not to be particularly significant, with the square cell model 

performing only slightly worse than Equation (20). 1 


1 We also note in passing that sheared environments favor formation of salt sheets aligned in the 
direction of the background current (Linden, 1974; Kimura and Smyth, 2007), in which case salt finger 
dynamics become similar to the 2D case discussed in Sec. 3a. It is perhaps ironic that since large-scale 
shears are ubiquitous in the ocean, salt fingers may be better represented by two-dimensional than by three- 
dimensional simulations. 
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The calculation of secondary growth rates for the square salt fingers Equation (20) 
proceeds, as previously, by the method based on the Floquet theory. The governing 
equations, Equation (2), are linearized with respect to the basic state Equation (18); 
which is assumed to be maintained in a quasi-steady state, and the solution is sought in 
the following form: 




C j 1 ^ 

n x’ n y 

S' 


S n„n y 

u 

.J 

N N 

u 

= Qxp(if x kx + if y ky + imz + At ) ^ ^ 

MX ’My 

V 

n x =-N n y =-N 

V 

n x’ n y 

W 


w 

yP'y 


^ , n y 


exp (in x kx + in y ky ), 


( 21 ) 


where (f x ,f y ) are the Floquet factors in x and y. Substituting Equation (21) in the 
linearized governing equations and collecting the individual Fourier components reduces 
the stability problem to matrix form Equation (13). Maximizing the growth rates with 
respect to (f x ,f y ,m), we predict the largest growth rates of secondary instabilities (m) as a 
function of (R p ,Pr,x,y4 T )- Finally, for any given C, we search for the amplitude A r = A* 
that satisfies the growth rate balance Equation (1). The vertical fluxes are then 
reconstructed as follows: 

F T =^T = 5 -4(\ + k 2 m „) 

' F -o 5 4(4+*L) 2 (22) 

= wS =---—— 

[ (V<,) 


and plotted in Figure 9 as a function of density ratio (R p ) for various values of C. 
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p 

Figure 9. Heat flux in the theoretical three-dimensional model for various values of C. 

While the patterns of 3D (Figure 9) and 2D (Figure 8) fluxes are qualitatively 
similar—in both cases fluxes rapidly reduce with R p and increase with C —we note that 
the 3D fluxes are consistently higher. The reason is attributed to the fact that the 2D 
basic state Equation (9) is more unstable than its 3D counterparts Equation (18) and 
Equation (20). Thus, the amplitude of salt fingers has to be significantly higher in 3D to 
produce secondary instabilities of the same strength as in 2D, which rationalizes the 
behavior noticed in the numerical simulations of Stern et ah, (2001). 
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IV. OCEANOGRAPHIC (HEAT-SALT) CASE: DNS ANALYSIS 


In this chapter, we present numerical simulations of oceanographic double- 
diffusive convection, warm saltier water over colder fresher water, with density ratio (R p ) 
values 1.1, 1.3, ..., 2.9 (Table 2). We have initiated these experiments by computer¬ 
generated small amplitude temperature and salinity distribution. Although double- 
diffusive interactions may lead to larger scale phenomena, such as thermohaline 
staircases, thermohaline intrusions, and internal waves (Traxler et ah, 2011), this 
experiment only focuses on 10 individual fingers (5 warm/salty and 5 cool/fresh 
initially). The vertical salinity (Fs) and temperature (Ft) fluxes were recorded as a 
function of time and their statistical averages were analyzed as a function of density ratio. 
Simulations of this nature have only recently become possible, as computer technology 
has exceeded the petaflop barrier. This has allowed us to explore other parameter regimes 
besides just heat and salt interactions as discussed in Chapter V. 

A. CALIBRATION OF THE THEORETICAL MODEL 

A limitation of the theoretical model in Chapter III is related to its inability to 
determine the constant C internally. Our theory assumes that C is comparable to unity 
but does not provide an exact value. In order to calibrate C, we now turn to DNS. For 
the heat/salt case (Pr=7,T=0.01) discussed in this section, the analysis is based on 10 two- 
dimensional and 10 three-dimensional simulations conducted for R p = 1.1,1.3,. ..,2.9. 

Figure 10 presents fluxes from a series of two-dimensional calculations, plotted 
along with the theoretical prediction from the foregoing growth rate balance theory 
(Chapter II) with C=4.3. In all simulations we used (N X ,N Z )=(256,512) grid points and 
the size of the computational domain (L X ,L Z )=(38.2,76.3) equivalent to 5 x 10 linearly 
fastest growing linger wavelengths ( d ). The variation of heat (Figure 10a) and salt 
(Figure 10b) fluxes with density ratio follows the theoretical pattern very closely. The 
ability of our model, containing only one adjustable parameter, to capture the details of 
the numerical simulations to such an extent lends credence to the proposed condition 

Equation (1) as the most relevant physical factor constraining growth of salt fingers. 
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Figure 10. Comparison of the theoretical two-dimensional model for C=4.3 (solid curve) 
with DNS (plus signs). Heat and salt fluxes are shown in a) and b ) respectively. 

Figure 11 presents the corresponding three-dimensional calculations. The 
resolution of three-dimensional DNS matched the two-dimensional runs: 
(N x ,N y ,N z )=(256,256,512) mesh was used to resolve 5x5x10 fastest growing finger 
wavelengths. The theoretical calculation for C=2.7 is superimposed on the numerical 
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data, revealing their mutual consistency—our theory predicts the right order of magnitude 
and the pattern of fluxes as a function of R p . The relevant value of the adjustable 
coefficient C=2.7 in 3D is comparable to—but differs from—the 2D value of C=4.3. The 
(limited) variation in C can be readily attributed to distinct equilibrium dynamics of 
three-dimensional and two-dimensional salt fingers in the unbounded model (Radko and 
Stern 1999). 




Figure 11. The same as in Figure 10 but for the three-dimensional model with 0=2.7. 
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B. 


OBSERVED PATTERNS 


T(x,y) S(x,y) 



Figure 12. Horizontal distribution of temperature and salinity for fully equilibrated state. 

Oceanographic values of Pr = 7, x = 0.01, R p = 1.9. 

The horizontal and vertical structure of individual fingers can be clearly seen from 
Figures 12 and 13, respectfully. It is interesting to note that the temperature field is 
characterized by much larger scales than that of salinity. It is also evident that while both 
the temperature and salinity distributions are spatially correlated there is more variability 
within the salinity field. This latter feature is due to the fact that salinity diffuses on 
much smaller scales (100 times slower) than temperature. 
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Figure 13. Vertical distribution of temperature and salinity for the fully equilibrated state. 

Oceanographic values of Pr = 7, x = 0.01, R p = 1.9. 


31 



















c. 


MIXING CHARACTERISTICS 


Equilibrium Transport - Heat/Salt - Rhol .9 



Figure 14. Temperature and salinity fluxes as a function of time for the heat-salt case. Ft 
and F s values are calculated after equilibrium is reached. 

Temperature and salinity fluxes for a typical North Atlantic oceanographic case 
(Pr = 7, x = 0.01, R p = 1.9) are plotted in Figure 14. Note that fluxes in Figure 14 are 
given in dimensional units and the sign convention reflects downward transport of heat 
and salt (unlike Figure 5). Our model starts using random noise until equilibrium is 
reached and the data for oceanographic applications are taken from the fully equilibrated 
stage (lower panel of Figure 14). 

Using data similar to that in Figure 14 for a range of R p , we can calculate a mean 
of our fluxes after equilibrium is reached. This corresponds to typical conditions in the 
Atlantic thermocline for a range of density ratio values (R p ). These values are included in 
Table 2, along with the dimensional temperature and salinity diffusivities (Kjdim and 
Ksdim) which are of the same order of magnitude expected in the ocean. Also included in 

F 

Table 2 are Buoyancy Flux (Fs - Ft) and Flux Ratio (y = —) values for the heat salt case. 
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Heat/Salt 

256x256x512 

Rp 

f t 

F s 

Krdim(m A 2/s) 

Ksdim(m A 2/s) 

Buoyancy 

Flux 

Ratio 

l.i 

-253.3917 

-369.8733 

-3.471E-05 

-5.574E-05 

116.4816 

0.685077025 

1.3 

-114.1179 

-186.4594 

-1.563E-05 

-3.321E-05 

72.3415 

0.61202546 

1.5 

-70.2319 

-121.8825 

-9.622E-06 

-2.505E-05 

51.6506 

0.576226284 

1.7 

-49.5589 

-90.2265 

-6.790E-06 

-2.101E-05 

40.6676 

0.54927211 

1.9 

-40.8760 

-75.8794 

-5.600E-06 

-1.975E-05 

35.0034 

0.538696932 

2.1 

-35.0073 

-65.8659 

-4.796E-06 

-1.895E-05 

30.8586 

0.531493535 

2.3 

-28.9132 

-56.0535 

-3.961E-06 

-1.766E-05 

27.1403 

0.515814356 

2.5 

-27.3972 

-52.2506 

-3.753E-06 

-1.790E-05 

24.8534 

0.524342304 

2.7 

-23.0670 

-45.2340 

-3.160E-06 

-1.673E-05 

22.1670 

0.509948269 

2.9 

-22.7098 

-44.1900 

-3.11 IE-06 

-1.756E-05 

21.4802 

0.51391265 


Table 2. Heat - Salt temperatures (Ft) and salinity (Fs) fluxes for a range of density 
ratios, Rp. Also shown are the dimensionalized temperature and salinity 
diffusivities (K Tdim , Ksdim), Buoyancy Flux (Fs-F T ), and Flux Ratio (F T /Fs). 


While development of explicit parameterizations of double diffusion for large- 
scale ocean modeling is outside of our immediate goal, it could prove beneficial for more 
applied studies to note that the numerical data in Figure 11 can be accurately 
approximated by the following algebraic expressions: 


F s = ° s +b s , (a s ,b s ) = (135.7,-62.75) 

A> ~ 1 

r = a y exp (~b r R p ) + c y , ( a y ,b y ,c y ) = (2.709,2.513,0.5128) 
F t = yF s , 


(23) 


which are indicated in Figure 15 by the red curve. The structure of the expression for 
Fs(Rp) is suggested by theoretical arguments in Radko (2008), whereas the expression for 
the flux ratio y(R p ) represents a purely empirical exponential fit to the numerical data. 
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Figure 15. The best fit for Heat Salt (Pr=7; x=0.01) utilizing cftool in MATLAB to solve 

Equation (23). 
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V. GENERALIZATIONS 


In this section, we go beyond the oceanographically motivated heat-salt 
(Pr,x)=(7,0.01) case and test the generality of our conclusions by exploring a broader 
parameter range. The discussion in this section is based on the comparison of four 
representative series of three-dimensional numerical experiments with corresponding 
theoretical results: 

• Case 1: Low Prandtl number (Pr,x)=(0.1,0.01); 

• Case 2: High Prandtl number (Pr,x)=(100,0.01); 

• Case 3: High diffusivity ratio (Pr,x)=(7,0.1); 

• Case 4: High Prandtl number and High diffusivity ratio(Pr,x)=(1000,1/3); 

Unfortunately, DNS for more extreme parameter values (Pr—»co, Pr—>0, x—>0) 
become computationally prohibitive in 3D. However, our selection represents a 
substantial fraction of the numerically accessible range. In all cases, we performed a 
series of calculations with R p = 1.1,1.3,...,2.9 in computational domains corresponding to 
5 x 5 x 10 fastest growing finger wavelengths, recording the average equilibrium fluxes. 

Comparing Case 1 with our baseline heat-salt configuration allows us to examine 
specifics of the low Prandtl number regime, which is relevant for astrophysical 
applications. Likewise, differences between Case 2 and heat-salt provide some glimpse 
into the large Prandtl number regime, exemplified by double-diffusion in magma 
chambers. Case 3 differs from the heat-salt calculation by diffusivity ratio, which also 
varies dramatically between various applications—from x = 0.8 for heat/humidity in the 
atmosphere to x~T0' 6 in stellar interiors. Finally, Case 4 (Sugar Salt parameters) is a 
combination of an increase in Prandtl number and diffusivity ratio and behaves 
differently than the majority of runs throughout our experiments, but is referenced to 
show how these differences in parameters affect results. 
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Density Ratio 

R p = 1.9 

Pr 

X 

F t 

F s 

Heat-Salt 

7 

0.01 

-40.8760 

-75.8794 

Case 1 

0.1 

0.01 

-1.9140 

-4.6501 

Case 2 

100 

0.01 

-220.1848 

-364.8436 

Case 3 

7 

0.1 

-39.3424 

-60.2595 

Case 4 

1000 

1/3 

-769.4803 

-885.3976 


Table 3. Dependence of the equilibrium T-S fluxes on the molecular characteristics 
(Pr,x) in three-dimensional simulations with R p =1.9. 


A. OBSERVED PATTERNS 

We can gain a qualitative understanding of the effects of parameter changes by 
reviewing Figures 16 to 19. We note that for each of these figures, the equilibrium point 
has already been reached. 
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Case 1 


Case 2 



Pr=7;t=0.1 Pr=1000; x=l/3 

Figure 16. DNS model of horizontal temperature (highest rate of diffusion) in (x, y) 

coordinates for Cases 1 through 4. 

Taking a horizontal slice of the temperature field (Figure 16) allows us to further 
explore the similarities and changes in fingers caused by differences in background 
parameters. It can be seen that there exists decreasing dependence between Prandtl 
number and horizontal finger scales. 
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Case 1 


Case 2 




Pi=7;t=0.1 



I 


I 


Figure 17. DNS model of horizontal salinity (least rate of diffusion) in (x,y) coordinates 

for Cases 1 through 4. 


At first glance, a review of the horizontal sections of salinity fields (Figure 17) 
does not show as clear cut of a correlation as that of temperature fields. By breaking the 
analysis down into groups, where we hold some parameters constant, we begin to see a 
pattern emerge. For example, from comparison of heat-salt (Figure 12) with cases 1 and 
2, where x = 0.01 and the Prandtl number changes from 7, 0.1, to 100, it can be seen that 
indeed diameter size increases with a decrease in Prandtl number. We can also examine 
cases where the Prandtl number remains constant and the diffusivity ratio varies, by 
comparing heat-salt and Case 3 (where x increases from 0.01 to 0.1). It is readily 
observed that increased x results in increased finger scales. 
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Figure 18. DNS model of vertical temperature (fastest rate of diffusion) in (x, z) 

coordinates. 
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In addition to what has been observed in the temperature field (horizontal) section 
example (Figure 16), we can now examine finger length scales and their dependencies 
from the vertical sections of the temperature fields (Figure 18). It is clear that thinner, 
longer fingers occur for higher Prandtl numbers. This graphic shows that a decrease in 
Prandtl number and an increase in the diffusivity ratio both cause a less defined finger 
formation. 


40 



Salinity (x.y) 
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Figure 19. DNS model of vertical salinity (slowest rate of diffusion) in (x, z) coordinates. 
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Visual inspection of the vertical salinity sections in Figure 19 supports the 
foregoing inferences. Comparing heat-salt with Cases 2 and 4 implies that an increase in 
Prandtl number results in thinner, longer individual fingers. Magma chamber fingers 
(with large Prandtl values) have been observed in laboratory experiments as being slim or 
slender (Tait 1989) and is verified by these results. A more objective measure of the 
previously observed patterns of salt fingers is given by the aspect ratio: 


A r = 


<T 2 +T(. +T 7 2 > 
3<T z 2 > 


(24) 


This quantity reduces to one in the limit of isotropic fingers and increases as fingers 
become more elongated. Comparison of heat salt and Case 3 (one order of magnitude 
increase in diffusivity ratio) shows only a slight increase in the aspect ratio (Ar). A 
summary of the dependencies of A R is shown in Figure 20 and in Table 4. 



Figure 20. Aspect Ratio, Equation (24), for R p = 1.9 for Heat Salt and Cases 1 to 4. 
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Aspect Ratio at R 

p=1.9 

Heat Salt 

Case 1 

Case 2 

Case 3 

Case 4 

1.4928 

1.4349 

6.5011 

1.7506 

12.5670 


Table 4. Aspect Ratio of Heat - Salt and Cases 1 through 4. 


For the same density ratio, R p = 1.9, (Figure 20 and Table 4) the aspect ratio 
increases with an increase in Prandtl number. As x increases, the aspect ratio increases 
slightly as reflected by the differences between the heat salt case and Case 3. 



Figure 21. Case 1 Aspect Ratio as a function of Density Ratio. 


Examination of Case 1 (Figure 21 and Table 5), in which Pr and x are held 
constant throughout a wide range of R p , reveals that an increase in density ratio is 
accompanied by a mono tonic increase in aspect ratio. 
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Case 1 

R P A r 

1.1 

1.2261 

1.3 

1.3271 

1.5 

1.3453 

1.7 

1.4101 

1.9 

1.4349 

2.1 

1.4796 

2.3 

1.4489 

2.5 

1.5791 

2.7 

1.6038 

2.9 

1.6564 


Table 5. Aspect Ratio as a function of density ratio values in Case 1. 
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B. 


TRANSPORT CHARACTERISTICS 


Equilibrium Transport - Case 1 - Rhol.9 



Equilibrium Transport- Case 4 - Rhol.9 



Figure 22. Temperature and salinity fluxes as a function of time for cases 1 through 4. 
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Case #1 

Rp 

Ft 

F s 

Buoyancy 

Ratio 

1.1 

-6.4812 

-14.4047 

7.9235 

0.449936 

1.3 

-4.1979 

-9.8905 

5.6926 

0.424438 

1.5 

-3.1962 

-7.4937 

4.2975 

0.426518 

1.7 

-2.581 

-6.0696 

3.4886 

0.425234 

1.9 

-1.914 

-4.6501 

2.7361 

0.411604 

2.1 

-1.6215 

-3.8882 

2.2667 

0.417031 

2.3 

-1.4437 

-3.4342 

1.9905 

0.420389 

2.5 

-1.2994 

-3.0476 

1.7482 

0.426368 

2.7 

-1.0685 

-2.5368 

1.4683 

0.4212 

2.9 

-0.9298 

-2.1956 

1.2658 

0.423483 


Case #2 

Rp 

Ft 

F s 

Buoyancy 

Ratio 

1.1 

-1105.08 

-1478.96 

373.8807 

0.7472 

1.3 

-468.524 

-705.2 

236.6764 

0.664384 

1.5 

-315.245 

-499.912 

184.6674 

0.6306 

1.7 

-242.235 

-398.488 

156.2523 

0.607887 

1.9 

-220.185 

-364.844 

144.6589 

0.603505 

2.1 

-196.214 

-328.582 

132.3677 

0.597154 

2.3 

-169.363 

-282.033 

112.6701 

0.600507 

2.5 

-168.507 

-286.134 

117.6272 

0.588908 

2.7 

-144.364 

-245.818 

101.4536 

0.587281 

2.9 

-149.324 

-250.356 

101.0315 

0.596448 


Case #4 

Rp 

Ft 

Fs 

Buoyancy 

Ratio 

1.1 

-3009.76 

-3407.43 

397.6623 

0.883295 

1.3 

-1112.64 

-1321.63 

208.9909 

0.841869 

1.5 

-873.25 

-1047.29 

174.0426 

0.833817 

1.7 

-870.275 

-1032.55 

162.271 

0.842844 

1.9 

-769.48 

-885.398 

115.9173 

0.869079 

2.1 

-868.552 

-983.512 

114.9605 

0.883112 

2.3 

-650.989 

-720.552 

69.5631 

0.903459 

2.5 

-578.38 

-626.742 

48.3624 

0.922835 

2.7 

-328.609 

-344.985 

16.3752 

0.952534 

2.9 

-14.0786 

-14.4201 

0.3415 

0.976318 


Case #3 

Rp 

Ft 

Fs 

Buoyancy 

Ratio 

1.1 

-234.914 

-324.005 

89.0909 

0.725032 

1.3 

-104.725 

-150.178 

45.453 

0.697339 

1.5 

-66.5649 

-98.8597 

32.2948 

0.673327 

1.7 

-49.8729 

-75.2535 

25.3806 

0.662732 

1.9 

-39.3424 

-60.2595 

20.9171 

0.652883 

2.1 

-33.3554 

-51.4704 

18.115 

0.64805 

2.3 

-26.356 

-41.2481 

14.8921 

0.638963 

2.5 

-12.2223 

-20.0053 

7.783 

0.610953 

2.7 

-15.3988 

-25.0506 

9.6518 

0.614708 

2.9 

-12.0376 

-19.733 

7.6954 

0.610024 


Table 6. Tables of DNS results for all cases including average temperatures (F T ) 
and salinity (Fs) fluxes, buoyancy (Fs - F T ) and flux ratio (F T /Fs). 


Typical evolutionary patterns of temperature and salinity fluxes (not necessarily 
representing temperature and salinity) are presented in Figure 22, for R p = 1.9. The 
summary of all calculations in Table 6 indicate that fluxes are most strongly affected by 
variation in the density ratio (R p ). Our limited numerical exploration of the effects of the 
molecular parameters suggests that, for any given R p , fluxes increase with the Prandtl 
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number (Pr) and very mildly decrease with the diffusivity ratio (x)—see Tables 3 and 6. 
These dependencies are consistent with the prediction based on the growth rate balance 
model. 

The purpose of the following calculation is to detennine whether the empirical 
parameterization (Equation (23)) developed for the oceanographic heat salt case applies 
equally well to other configurations. This parameterization was applied to each case, 
using the cftool in MATLAB. Results are summarized in Table 7 (with quality of the fit 
indicated in the last column) and Figure 23. 


Empirical Parameterization 


a s 

b s 

a Y 

b v 

c v 

Goodness 
Of Fit 

Heat Salt 

135.7 

-62.75 

2.708 

2.513 

0.5128 

99.7% 

Case 1 

5.113 

-0.7788 

427.7 

8.738 

0.4212 

95.9% 

Case 2 

503.9 

-160.6 

8.649 

3.663 

0.5928 

99.0% 

Case 3 

124.6 

-72.75 

0.359 

0.8303 

0.5772 

99.8% 

Case 4 

1189 

-482.2 

N/A 

N/A 

N/A 

93.5% 


Table 7. Empirical parameterization of the temperature and salinity fluxes. 
Parameterization of the flux ratio for Case 4 has not been attempted due to the 
nonmonotonic pattern of y(R p ), which is inconsistent with the assumed expression, 

Equation (23). 


Results in all cases are encouraging, suggesting the universal character of the 
assumed parameterization similar to Equation (23), as quantified by the goodness of fit 
for all simulations. Examining Table 7 shows the Prandtl number value has a significant 
effect on the magnitude of the a s and b s coefficients: as Prandtl number increases so do 
the magnitude of these values. 
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Figure 23. Best fit plots for Table 7, Cases 1 through 4, respectively. Compare to Fleat 

Salt Best Fit (Figure 15). 
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Buoyancy Flux 
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Figure 24. Buoyancy flux (Fs - Ft) is plotted showing that as R p increases buoyancy 
decreases. This is an indirectly proportional relationship. 


One of the key characteristics of fully developed salt fingering is related to its 
ability to transport buoyancy in the vertical direction. This characteristic is thought to be 
involved in the generation and maintenance of thermohaline staircases (Merryfield, 
2000). The buoyancy flux values of the cases presented in this study (Fs - F T ) are 
examined in Figure 24, which reveals that an increase in R p results in a decrease in the 
buoyancy. 
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Buoyancy Flux 



Figure 25. Buoyancy flux for each case at R p = 1.9. 


By taking a closer look at all cases at a fixed density ratio, R p = 1.9, (Figure 25 
and Table 8) we can see the two largest Prandtl number values in cases 2 and 4 (Pr =100 
and Pr = 1000 respectively) have the largest buoyancy values. In comparing the heat salt 
and Case 3 in which only the diffusivity ratio changes (x = 0.01 and x = 0.1, respectively) 
there is a slightly smaller buoyancy value for Case 3. Finally, in comparison of heat salt 
and Case 1, with a decrease in Prandtl number, the buoyancy decreases. 


Buoyancy flux at R n =1.9 

Heat Salt 

Case 1 

Case 2 

Case 3 

Case 4 

35.0034 

2.7361 

144.6589 

20.9171 

115.9173 


Table 8. Buoyancy flux at R p = 1.9 of all cases. 
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Figure 26. Flux ratio in all case studies. 


Another critical parameter, also implicated in the formation of thermohaline 


staircases, is given by the flux ratio 


Y = - 


. The flux ratio generally decreases with an 


L s y 


increase in R p as seen in Figure 26. However, the two cases (3 and 4) with larger 
diffusivity ratio (t) show an increase in the flux ratio for higher values of density ratio 
(R p ). Figure 27 and Table 9 represent the flux ratio for each case at a fixed value R p = 
1.9, which also indicate that the largest values of flux ratios occur for the largest value of 
diffusivity ratio. 
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Flux Ratio 



Figure 27. Flux Ratio for all cases at R p = 1.9. 


Flux Ratio at R n =1.9 

Heat Salt 

Case 1 

Case 2 

Case 3 

Case 4 

0.538696 

0.4116 

0.603505 

0.65288 

0.8690788 


Table 9. Flux Ratio for all cases at R p = 1.9. 


C. COMPARISON WITH THE THEORETICAL MODEL 

The relevant values of the adjustable coefficient (Q vary only weakly between 
various cases considered in this study (see Figures 28 through 30). The lowest value of 
0=1.6 (obtained for Case 1) is not too different from that largest value of C=2.7 (heat-salt 
case) considering the range of governing parameters and fluxes. Thus, numerical 
evidence in this study supports the growth rate balance theory as a conceptual model of 
the equilibrium salt fingering, capable of predicting dependencies of the equilibrium 
temperature / salinity transports within a factor of two. 
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Figure 28. Case 1: Temperature (left panel) and salinity (right panel) fluxes as a function 
of density ratio (Pr= 0.1, x = 0.01). Solid curves indicate prediction of the 
theoretical three-dimensional model and DNS values are represented by plus signs. 



Figure 29. Case 2: Temperature (left panel) and salinity (right panel) fluxes as a function 
of density ratio (Pr= 100, t = 0.01). Solid curves indicate prediction of the 
theoretical three-dimensional model and DNS values are represented by plus signs. 
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Figure 30. Case 3: Temperature (left panel) and salinity (right panel) fluxes as a function 
of density ratio (Pr= 7, x = 0.1). Solid curves indicate prediction of the theoretical 
three-dimensional model and DNS values are represented by plus signs. 
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VI. DISCUSSION AND CONCLUSIONS 


Double-diffusion is an essential mixing process that may very well drive the 
larger scale oceanographic circulations and a more thorough study of this phenomenon 
may lead to higher accuracy in climate and ice melt models and global and regional 
oceanographic forecast models. Understanding this process requires a series of direct 
numerical simulations that would not have been possible without parallel processing. 
The processing time utilized in our study totals over 800,000 hours—350,000 hrs on 
Department of Defense (DoD) Systems and the remaining on NPS computers. In other 
words, it would take almost a full century to complete these calculations on a single 
processor. The parallelization of our de-aliased pseudo-spectral model allowed us to 
utilize a message passing interface to effectively resolve the complex calculations of 
fluxes, and therefore derive an understanding of how these micro-scale mixing processes 
occur in nature. Such an understanding will facilitate our future ability to model these 
interactions and their effects on large scale circulations. This too requires increased 
speeds of computer processors in order to process even more data within a timely, 
forecastable time period. 

Our model predicts the following. Double-diffusive transport of heat and salt 
rapidly intensifies with decreasing density ratio. Fluxes are less sensitive to molecular 
characteristics, mildly increasing with Prandtl number (Pr) and decreasing with 
diffusivity ratio (x). The aspect ratio of individual fingers monotonically increases with 
Prandtl number, with tall slender fingers realized for large Pr and irregular chaotic fingers 
for low Pr. Numerically deduced dependency of temperature and salinity fluxes on 
density ratios are adequately described by empirical Equation (23), with goodness of fit 
greater than 93% in all calculations. Flux ratio shows an increase at larger R p values for x 
>0.1 while constantly decreasing for x = 0.01. 

In addition to the insights based on the numerical modeling, this study presents a 
phenomenological theory for the vertical transport by salt fingers in unbounded 
temperature and salinity gradients, which is based on the growth rate balance Equation 
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(1). By utilizing DNS we were able to plot equilibrium fluxes of temperature and salinity 
as a function of time and explore their dependency on growth rates. This dependency led 
us to the discovery of the adjustable coefficient (C) that varies slightly for a wide range of 
parameters and fluxes. Theory assumes that the statistical equilibrium is reached when 
the growth rate of salt fingers becomes comparable to the growth rate of their secondary 
instabilities. This assumption can be rationalized as follows. If the growth rate of 
primary salt fingers exceeds that of secondary instabilities, then the primary instabilities 
are unlikely to be affected by perturbations and continue to grow. If, on the other hand, 
the secondary instabilities grow much faster than fingers, they would rapidly disrupt the 
primary modes, dramatically reducing their amplitude and transport characteristics. 
Thus, one is led to a conclusion that the approximate balance between the growth rates of 
primary and secondary instabilities is inevitable for statistical equilibrium and should be 
satisfied for a wide range of governing parameters. 

It should be understood that the growth rates, both primary (ki) and secondary 
(^ 2 ), can only be estimated in the order of magnitude sense. Calculation of A| is made for 
the fastest growing salt fingers - the vertically uniform elevator modes - whereas fully 
developed double-diffusive convection is generally dominated by irregular structures of 
finite vertical extent. Similar problems arise for secondary instabilities: iC is estimated 
for the steady and vertical background field, neither of which is an exact statement. In 
view of these difficulties, we only claim that the growth rates estimated by our theory 
should be comparable, but not necessarily equal. Thus, we consider a weak form of the 
growth rate balance in Equation (1). The nondimensional order one coefficient C in 
Equation (1) cannot be determined internally on the basis of theory and therefore it was 
calibrated using direct numerical simulations covering a wide range of governing 
parameters. For each (Pr,x), theory adequately describes the (rapidly decreasing) 
dependence of fluxes of diffusing components on density ratio. The calibrated 
coefficients C are fairly stable. In all three-dimensional simulations - simulations that 
cover three orders of magnitude in Pr and one order in x - the values of C are limited to a 


56 



relatively narrow range 1.6 < C < 2.7. The theoretical prediction for C = 2.2 can explain 
all the numerical fluxes in this paper within a factor of two. This can be readily 
rationalized using the growth rate balance Equation (1). 

Finally, it is important to emphasize the theoretical significance of our model. 
We believe that this study opens opportunities for purely analytical explorations that are 
focused on the growth rate balance Equation (1), rather than on the more complicated 
original Navier-Stokes system. Our optimism is partially based on the pronounced low- 
order dynamics of our model. Thus, it is likely that the general character of our solutions 
could be captured by an analytical low-dimensional model. 
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VII. RECOMMENDATIONS FOR FUTURE WORK 


Imagine our children and grandchildren planning experiments on Jupiter or other 
celestial bodies. The work we have done here could play a pivotal role in how they 
design their sensors and plan their missions. As we have no easy way to achieve direct 
measurements of these planetary systems, reliable data is often scarce and the spatial 
location of the elements that need to be measured must be pinpointed to the most 
scientifically relevant areas. Environmental characteristics of any future planetary 
system, need to be pinpointed to allow the maximum effectiveness of these future 
missions. 

Even closer to home imagine forecast models, atmospheric and oceanographic, 
with successful output out to 30 days and beyond. Use of a simple, computationally 
compact, algorithm for parameterization of micro-scale mixing in general and salt fingers 
in particular could be implemented in future forecast models would save both processing 
time and money. Further exploitation of this microscale feature and its impact on large 
scale circulations can greatly improve climate and operational oceanographic models. Ice 
melt and climate models could also gain accuracy with the implementation of double 
diffusive convection. As computer technology improves and processing speeds increase 
this less computationally taxing algorithm can easily be input into these forecast models. 

Successfully modeling double diffusive convection could impact models on 
magma chamber flows and may even prove vital in volcanic eruptions and possible 
earthquake predictions. With successful earthquake and volcanic eruption predictions we 
can also predict tsunamis and have more warning time in order to prepare for these often 
catastrophic events. The events in the past few years show that these advanced warning 
systems are needed to prevent the mass loss of life with the “Ring of Fire” surrounding 
the Pacific Ocean. 
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